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ABSTRACT 

We analyse the properties of substructures within dark matter haloes (subhalos) using 
a set of high-resolution numerical simulations of the formation of structure in a ACDM 
Universe. Our simulation set includes 11 high-resolution simulations of massive clus- 
ters as well as a region of mean density, allowing us to study the spatial and mass 
distribution of substructures down to a mass resolution limit of 1O 9 /i _1 M0. We also 
investigate how the properties of substructures vary as a function of the mass of the 
'parent' halo in which they are located. We find that the substructure mass function 
depends at most weakly on the mass of the parent halo and is well described by a 
power-law. The radial number density profiles of substructures arc steeper in low mass 
haloes than in high mass haloes. More massive substructures tend to avoid the centres 
of haloes and are preferentially located in the external regions of their parent haloes. 
We also study the mass accretion and merging histories of substructures, which we 
find to be largely independent of environment. We find that a significant fraction of 
the substructures residing in clusters at the present day were accreted at redshifts 
z < 1. This implies that a significant fraction of present-day 'passive' cluster galaxies 
were still outside the cluster progenitor and were more active at z ~ 1. 
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1 INTRODUCTION 

The formation and evolution of structure in the Universe 
is a topic of fundamental interest. In the last decades, the 
Cold Dark Matter (CDM) model has been extensively stud- 
ied and has had considerable success in reproducing obser- 
vational results, both on galactic and on cluster scales. In 
fact, the CDM model with the 'concordance' set of cosmo- 
logical parameters (ACDM) has been so successful that it 
can now be considered a standard paradigm for the forma- 
tion of structure in the Universe. According to this model, 
the dominant force that drives structure formation is grav- 
ity, and large systems like galaxy clusters are formed via 
hierarchical merging of smaller structures. 

Numerical simulations of gravitational clustering of 
dark matter are an indispensable tool for investigating the 
non-linear growth of structures in its full geometrical com- 
plexity. Until recently, dissipationless simulations suffered 
from the so-called overmerging problem, i.e. substructures 
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disrupt very quickly within dense environments (Katz & 
White 1993). However, both analytic work (Moore, Katz 
& Lake 1996) and high resolution simulations (Tormen, 
Bouchet & White 1997; Ghigna et al. 1998; Klypin et al. 
1999b; Ghigna et al. 2000) have demonstrated that the cores 
of dark matter haloes that fall into a cluster can actually 
survive as self-gravitating objects orbiting in the smooth 
dark matter background of the cluster, provided high enough 
force and mass resolution are used. Recent high-resolution 
simulations (Ghigna et al. 2000; Springel et al. 2001b) have 
also shown that the abundance of these substructures is in 
agreement with the observed abundance of galaxies in clus- 
ters, suggesting a natural one-to-one identification of lumi- 
nous cluster galaxies and dark matter substructures. 

Another interesting claim is that the shape of the sub- 
structure mass function is independent of the mass of the 
parent halo (Moore et al. 1999b). It is not obvious that this 
should be the case, because in CDM cosmologies, the initial 
conditions do depend on scale and galaxies form several bil- 
lion years before clusters. As Moore et al. (1999b) pointed 
out, the logarithmic slope of the power spectrum asymptot- 
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ically approaches —3 on small scales, so clumps of widely 
different (yet sufficiently small) mass tend to collapse simul- 
taneously, and as a result the timescale between the collapse 
of the first substructures and their incorporation into larger 
haloes becomes shorter. One might therefore expect that 
substructures were more easily disrupted in low mass haloes, 
producing a substructure mass function which depends on 
the mass of the parent halo. 

Observationally, the predicted abundance of substruc- 
tures in clusters is one of the major successes of the CDM 
model (Springel et al. 200fb), but on galactic scales, it 
appears that simulations predict more substructures than 
are visible by almost two orders of magnitude (Kauffmann, 
White & Guiderdoni 1993; Moore et al. 1999a; Klypin et al. 
1999b; Tasitsiomi 2002). This is commonly referred to as the 
'dwarf galaxy crisis' of CDM. There have been suggestions 
that the solution to this problem lies in processes such as 
heating by a photo-ionising background that suppresses star 
formation in small haloes at early times (Efstathiou 1992; 
Bullock, Kravtsov & Weinberg 2000; Somerville 2002; Ben- 
son et al. 2002). Alternatively, it has been suggested that 
the nature of dark matter may be different than assumed 
in the canonical ACDM model, for example by being warm 
or self -interacting, both of which could selectively eliminate 
small-scale structure. However, self-interactions appear to 
be relatively ineffective in reducing the number of subha- 
los, unless the assumed cross-section is unreasonably large 
(Colin et al. 2002). 

A less drastic resolution was suggested by Stoehr et 
al. (2002) who noted that it might be possible to identify 
the observed Galactic satellites with the few most massive 
subhalos and that the rest contain no stars. Direct evidence 
for the large population of dark satellites predicted by CDM 
models comes from the anomalous flux ratios of multiply 
imaged quasars (Mao & Schneider 1998; Chiba 2002; Dalai 
& Kochanck 2002). 

So far, a detailed numerical analysis of substructures 
has only been carried out in high resolution re-simulations 
of a few individual haloes (Moore et al. 1999b; Ghigna et 
al. 2000; Springel et al. 2001b). In this paper, we carry out 
a systematic analysis of substructures as a function of the 
mass of the parent halo and as a function of environment. 
We study the mass functions of subhalos, their radial dis- 
tributions and their merging and mass accretion histories. 
These quantities are of fundamental interest for galaxy for- 
mation, because dark matter haloes and substructures rep- 
resent the birth places of luminous galaxies. Their accretion 
and merging histories regulate the rate at which baryons can 
cool, determining in this way the rate at which stars form 
in galaxies as a function of cosmic time. 

The layout of the paper is as follows: in Section 2 we 
describe the simulations that we use; in Section 3 we give a 
brief description of the algorithm used to find substructures 
in haloes; in Section 4 we present the subhalo mass function, 
in Section 5 we analyse the mass distribution of the largest 
substructures; in Section 6 we study the radial distribution 
of substructures; in Section 7 we discuss the merging and 
mass accretion histories of substructures, both as a function 
of mass and as a function of environment. A summary and a 
discussion of the results obtained are presented in Section 8. 



2 CLUSTER SIMULATIONS 

In this study, we use collisionless simulations of clusters gen- 
erated using the 'zoom' technique (Tormen et al. 1997; Katz 
& White 1993). First, a cosmological simulation of a large 
region is used to select a suitable target cluster. The parti- 
cles in the final cluster and its surroundings are then traced 
back to their initial Lagrangian region and are replaced by a 
larger number of lower mass particles. These are perturbed 
using the same fluctuation distribution as in the parent sim- 
ulation, but now extended to smaller scales to account for 
the increase in resolution. This resampling of the initial con- 
ditions of the Lagrangian region of the cluster thus allows a 
localised increase in mass and force resolution. Outside the 
high-resolution region, particles of variable mass, increasing 
with distance, are used so that the computational effort is 
concentrated on the cluster of interest, while still maintain- 
ing a faithful representation of the large-scale density and 
velocity field of the parent simulation. 

In this paper, we study a set of 11 high-resolution re- 
simulations of galaxy clusters (5 of mass 1O 14 /i _1 M and 
6 of mass 10 15 ft _1 Mg), and a high resolution re-simulation 
of a 'typical' region of the Universe. The simulations were 
carried out with the parallel tree-code GADGET (Springel, 
Yoshida & White 2001). 

One of our massive clusters was taken from the 'S-series' 
studied by Springel et al. (2001b), where the parent simu- 
lation employed was the GIF-ACDM model carried out by 
the Virgo Consortium (Kauffmann et al. 1999). This parent 
simulation followed 256 3 particles of mass 1.4 x 10 10 ft _1 MQ 
within a comoving box of size 141.3 /i _1 Mpc on a side. The 
other cluster re-simulations and the simulation of the field 
region were selected from the VLS simulation carried out by 
the Virgo Consortium (Jenkins et al. 2001; Yoshida, Sheth & 
Diaferio 2001). The simulation was performed using a par- 
allel P3M code (Macfarland et al. 1998) and followed 512 3 
particles with a particle mass of 7x 10 10 hT 1 Mq in a comov- 
ing box of size 479 fe _1 Mpc on a side. In all cases, the par- 
ent simulation and the re-simulations were characterised by 
the following cosmological parameters: Qo — 0.3, Qa = 0.7, 
spectral shape F = 0.21, h — 0.7 (we adopt the convention 
Ho — 100/ikms _1 Mpc _1 ) and normalisation as = 0.9. 

In Table 1, we summarise some important numerical 
parameters of the simulations used. We will refer to our 
five high-mass clusters of mass resolution 2 x 10 9 /i^'Mq as 
type 'Bl', and to the low-mass clusters as type 'B2'. These 
simulations were carried out by Barbara Lanzoni as part of 
her PhD thesis and were previously used in Lanzoni et al. 
(2003). The 'S2' simulation is taken from the 'S-series' of 
Springel et al. (2001b). The field region 'M3' was adopted 
from the 'M-series' studied by Stoehr (2003). 



3 IDENTIFICATION OF DARK MATTER 
SUBSTRUCTURES 

The identification of substructures in dark matter haloes 
is a difficult technical problem and many different 
algorithms have been developed to accomplish this 
task, for example the hierarchical friends-of-friends al- 
gorithm (HFOF) (Gottlober, Klypin & Kravtsov 1998), 
the bound density maximum algorithm (BDM) (Klypin 
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Table 1. Numerical parameters for the simulations used. All the simulations were carried out assuming a ACDM cosmology with 
cosmological parameters Ho = 0.3, SIa = 0.7, T = 0.21, erg = 0.9, and h = 0.7. In the table, we give the particle mass m p in the high 
resolution region, the starting redshift z s tart i the gravitational softening e in the high-resolution region and the number of simulations in 
each group N. 



Name Description m p [h 1 Mq] z s tart £ [h x kpc] N 



Bl 


10 15 /i _1 M clusters 


2 X 10 9 


60 


5.0 


5 


B2 


10 14 h~ 1 M Q clusters 


2 X 10 9 


60 


5.0 


5 


S2 


10 15 ^- 1 M cluster 


1.36 x 10 9 


50 


3.0 


1 


M3 


field simulation 


1.7 x 10 8 


120 


1.4 


1 



ct al. 1999a), and SKID (see http ref: http://www- 
hpcc.astro.washington.edu/tools). Each of these algorithms 
has its own advantages and weaknesses, so that arguably 
none of them is completely satisfactory yet. In this work, 
we use the algorithm SUBFIND proposed by Springel et al. 
(2001b), which combines ideas used in other group finding 
techniques with a topological approach for finding substruc- 
ture candidates. SUBFIND can handle haloes of arbitrary 
shape, does not require an iterative procedure for finding 
subhalo candidates, and is capable of detecting arbitrary lev- 
els of 'subhalos within subhalos'. In this section, we briefly 
summarise how the method works. 

In a first step, a standard friends-of-friends (FOF) algo- 
rithm is used to identify virialized parent haloes. The FOF 
algorithm links together all particle pairs with separation 
less than a linking length b. We adopt the standard value 
b = 0.2 in units of the mean particle separation, which se- 
lects groups of particles with overdensities close to the value 
predicted by the spherical collapse model for the virialized 
regions of haloes. The next step is to compute an estimate 
of the local density at the position of each particle in the 
group. To this end, we employ an adaptive kernel interpo- 
lation method similar to the one used in smoothed particle 
hydrodynamics. In the resulting density field, we define as 
substructure candidates locally overdense regions which are 
enclosed by isodensity contours that traverse a saddle point. 
Our method for finding these regions can be visualised as 
follows: we reconstruct the density field by considering par- 
ticles in order of decreasing density, thus working our way 
from high to low density. This corresponds to gradually low- 
ering a global threshold in the density field sampled by the 
simulation particles. Isolated overdense regions grow slowly 
in size during this process. When two such separate regions 
coalesce to form a single region, their density contours join 
at a saddle point. Each time such an event occurs, we have 
found two substructure candidates. 

After the regions containing substructure candidates 
have been identified, we apply an unbinding procedure where 
we iteratively reject all particles with positive total energy in 
order to eliminate 'background' particles that do not belong 
to the subhalo. For the purposes of this study, we consider 
all substructures that survive this unbinding procedure, and 
still have at least 10 self-bound particles, to be genuine sub- 
halos. 

In summary, the algorithm SUBFIND decomposes a 
given particle group into a set of disjoint and self-bound 
substructures, each of which is identified as a locally over- 
dense region in the density field of the background halo. 
Note that SUBFIND classifies all the particles inside a FOF 



group either as belonging to a bound substructure or as be- 
ing unbound. The self-bound part of the FOF background 
halo itself will then also appear in the substructure list. We 
will exclude it when referring to subhalos or substructures 
in the following analysis. 



4 THE SUBHALO MASS FUNCTION 

The sample of parent haloes used for studying the mass 
function analysis consists of 6, 5, 34 and 100 haloes in the 
mass ranges 8.68 x 10 14 -1.79 x 10 15 (from simulations Bl 
and S2), 6.99 x 10 13 -1.27x 10 14 (from simulations B2), 7.0 x 
10 12 -2.0 x 10 13 (from simulation M3) and 7.0 x 10 n -2.0 x 
10 12 h~ 1 M e (from simulation M3). 

The resulting subhalo mass functions are shown in 
Fig. 1. In the first four panels, we plot the differential mass 
functions for parent haloes of different mass. The histograms 
are computed by stacking all the haloes in the given range of 
mass and the error bars represent Poisson errors. The solid 
line in each of the panels is a power-law fit to the measured 
differential mass function; the fit is performed using the least 
absolute deviation method over the range of mass shown by 
the line. In all the cases the slope of this unrestricted fit is 
close to —1 (it is equal to —0.98 for the top left panel, —0.97 
for the top right panel, —1.11 for the middle left panel and 
— 1.13 for the middle right panel). However, we note that 
the lowest mass bins, which have the smallest statistical er- 
rors, are best fit with a slightly shallower slope: if we restrict 
the fit to the 4 lowest mass bins the slope is —0.94 for the 
top left panel and —0.85 for the middle left panel. These 
are closer to the value —0.8, measured by Helmi, White & 
Springel (2002) for a single cluster simulation of extremely 
high- resolution. Also note that a slope shallower than —1 at 
the low-mass end implies that the integrated mass in sub- 
structures remains bounded and is dominated by the most 
massive subhalos. It is likely that our subhalo mass func- 
tions are steepened somewhat by a cut-off in abundance for 
very massive substructures. 

The bottom left panel of Fig. 1 shows the cumulative 
mass function for all the haloes used in the sample. To com- 
pare the different subhalo mass functions, we have rescaled 
the subhalo mass by dividing by the virial mass of the parent 
halo. Each line represents the average cumulative mass func- 
tion over all the haloes in each mass bin. Note that in this 
paper we define the 'virial mass' as M200, the mass within 
a sphere of density 200 times the critical value at redshift 
zero. The lines end at different places because of the differ- 
ing mass resolution of the simulations (see Table 1). We find 
that all four cumulative mass functions agree within the star 
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Figure 2. The cumulative mass function of subhalos in units 
of rescaled subhalo mass for haloes with mass ~ 10 14 /i _1 Mq. 
The solid line corresponds to the average of haloes from simula- 
tions B2 and the dashed line to haloes from simulation M3. The 
vertical dotted line shows the resolution limit corresponding to 
simulations B2. 



tistical errors. Finally, in the bottom right panel, we show 
the differential subhalo mass functions in units of rescaled 
mass. 

We note that the 'universality' of the subhalo mass func- 
tion seen here appears to be quite robust with respect to 
numerical resolution. In Fig. 2 we compare the average cu- 
mulative mass functions for haloes with mass ~ 10 14 ft _1 M 
from simulations B2 and M3. We here averaged 5 haloes 
for simulations B2, and 4 for simulation M3, to reduce the 
object-to-object scatter that is unavoidable for small num- 
bers of subhalos. Despite an order of magnitude difference in 
numerical resolution, the agreement between the simulations 
is good. We are able to resolve the right number of objects 
in the low-resolution simulation above its resolution limit 
(shown as a vertical dotted line in the figure). A similar re- 
sult was obtained by Springel et al. (2001b, see their Fig. 5) 
using a set of 4 re-simulations of the same cluster with sys- 
tematically increasing resolution, thereby allowing a direct 
study of numerical convergence. This showed in particular 
that the S2 simulation used here has well converged to the 
properties of a much higher resolution simulation above its 
own resolution limit, as used here. Further support for our 
results was also found by Stoehr et al. (2003, see their Fig. 3). 
They compared the S-series simulations from Springel et al. 
(2001b) with an extremely well resolved re-simulation of a 
Milky-Way sized halo. This latter simulation used an up- 
dated version of the simulation code and more conservative 
integration parameters than used here (following Power et 
al. 2003), suggesting that the subhalo mass function is a rel- 
atively robust quantity and that the simulations we discuss 
here are adequate for our purposes. 

As a further check of the robustness of our results we 



also compare the internal structure of subhalos drawn from 
our different simulations. Fig. 3 shows the correlation be- 
tween the substructure mass and the third power of the 
maximum circular velocity, Vmax, measured directly from 
the circular velocity curve of each subhalo. Different sym- 
bols are used for substructures drawn from different simula- 
tions. Note that for the range of masses shown in the plot, 
substructures drawn from simulation M3 contain at least 60 
particles. While the scatter is clearly large for haloes with 
such a low number of particles, the good general agreement 
between the runs suggests that the smallest substructures in 
our lower resolution simulations have an internal structure 
that is still reliably resolved, at least in a statistical fashion. 

Our results confirm the conclusion drawn by Moore et 
al. (1999b): the mass function of substructures appears to be 
almost independent of the mass of the parent halo. While 
our results are consistent with such a 'scale-free' subhalo 
mass function, the halo-to-halo scatter in our simulation set 
is quite large, preventing us from putting tight constraints 
on the accuracy with which the 'scale-free' subhalo mass 
function is preserved when haloes of different mass are con- 
sidered. And hence there is still room for weak trends with 
mass. A clear detection of these would require simulations 
with larger dynamic range, and larger samples of simulated 
haloes for each mass bin. 

As we discuss in more detail in Sec. 7.2, our findings 
suggest that the destruction of satellites due to the physical 
processes of dynamical friction and tidal stripping on one 
hand, and the accretion of new satellites on the other hand, 
cancel out in such a way that the subhalo mass function does 
not depend or at best very weakly depends on the mass of 
the parent halo. The reason for the invariance of the subhalo 
mass function may lie in the physical nature of this dynam- 
ical balance, which may be insensitive to the slightly broken 
scale-invariance of dark haloes themselves. This shows up 
as a mass-dependence of halo concentrations, for example. 
Some fully analytic models for the subhalo abundance have 
been developed (e.g. Sheth 2003), but they are presently not 
able to account for mass-loss and dynamical friction self- 
consistently, and so provide little guidance to answer this 
interesting question. A full understanding of the apparent 
'conspiracy' that establishes an almost mass-invariant sub- 
halo mass function will therefore require further simulations. 



5 THE MOST MASSIVE SUBSTRUCTURES 

In this section we investigate whether the properties of the 
largest substructures depend on the mass of the parent halo. 
This is interesting since the largest substructures mark the 
sites where one expects to find the brightest galaxies. 

In the following, Mi refers to the mass of the most 
massive subhalo and M2 to the mass of the second most 
massive subhalo within the virial radius of a given object 
of virial mass M2oo- Note that we have excluded from our 
analysis the subhalo associated with the FOF group itself. 
In a semi-analytic scheme, it is this 'subhalo' that would 
host the brightest cluster galaxy (BCG). In Sec. 7.2 we will 
show that, once accreted onto a massive halo, substructures 
suffer significant stripping, an effect that is more important 
for substructures accreted at higher redshift. It is then likely 
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Figure 1. Panel (a) — (d): the differential mass function of subhalos residing in parent haloes of different mass. The solid line represents 
a power law fit to the mass function. The dashed line shown in panel (a) and (c) represents a fit restricted to the 4 lowest mass bins. 
The masses arc in units of /i- 1 M . The range of mass of the haloes used in each bin is indicated in each panel. Panel (e): the cumulative 
mass function of subhalos in units of rescaled subhalo mass. The solid line is for haloes with mass ~ 10 15 h~ 1 M.Q, the dashed line is for 
haloes with mass ~ 10 14 h~ 1 MQ, the dash-dotted line is for haloes with mass ~ 10 13 /i~ 1 Mq, and the long-dashed line is for haloes with 
mass ~ 10 12 /i -1 Mq. Panel (f): differential mass functions in units of rescaled subhalo mass; the different line styles are the same as in 
panel (c). 
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Figure 3. Substructure masses as a function of Vj5 ax for subha- 
los drawn from our different simulations. Small dots are used for 
subhalos that reside in haloes with mass ~ 10 15 h~ 1 Mq (from 
simulations Bl) and ~ 10 14 h^ 1 Mq (simulations B2) and filled 
circles for subhalos in haloes with mass ~ 10 13 h' 1 Mq (simula- 
tion M3). V m ax was determined as the maximum of the circular 
velocity curve of each subhalo. 



that the largest substructures we find within the virial radius 
at the present time were accreted at relatively low redshift. 

In Fig. 4, we plot M1/M200 as a function of M200 for 
434 haloes drawn from all the simulations listed in Table 1. 
This sample includes not only the central clusters in our 
re-simulations, but also the other haloes found in the high- 
resolution regions around the re-simulated objects down to a 
mass limit of 10 13 h" 1 Mq. We took care however to exclude 
contaminated haloes that contained low resolution particles. 
In simulation M3, we selected only haloes with a mass larger 
than 10 12 h~ 1 MQ. The small symbols in Fig. 4 indicate the 
value of M1/M200 measured for each individual halo, while 
the filled circles represent the median of the distribution. 
We have taken bins in M200 such that there are an equal 
number (143) of haloes in each bin, except for the last six 
points, which we treated as a separate bin, corresponding 
to the central cluster haloes in simulations Bl and S2. The 
error bars mark the 20 and 80 per cent percentiles of the 
distribution. 

The results in Fig. 4 suggest that M1/M200 depends 
very little on the mass of the parent halo. Interestingly, 
M1/M200 appears to exhibit less scatter for the most mas- 
sive haloes, but the number of simulated clusters we have in 
this high mass regime is rather small, so it is unclear whether 
this effect is statistically significant. 

The results in Fig. 4 imply that the median value of Mi 
increases in proportion to the mass of the cluster, suggest- 
ing that second ranked galaxies will be more luminous in 
more massive haloes. Note also that the mass of the largest 
substructure within the virial radius is typically only a few 
per cent of the virial mass. 



Figure 4. Ratio of the mass of the most massive substructure and 
the parent halo mass as a function of parent halo mass. The small 
symbols represent the values measured for each individual halo; 
the filled circles are the median in bins of halo mass chosen so 
that each of them contains the same number of points (143). The 
last six points, corresponding to the main haloes of simulations 
Bl and S2, are treated as a separate bin. The error bars mark the 
20 and 80 per cent percentiles of the distribution. 

In Fig. 5, we show the ratio M2/M1 as a function of 
the mass of the parent halo. Once again there is rather little 
dependence on M200 with a possible decrease in the scat- 
ter for more massive haloes. Note that the median value of 
M2/M1 is around 0.5. If the stellar masses of the second 
and third brightest galaxies in a cluster scale simply with 
the masses of their dark subhalos, they should have K-band 
luminosities that are equal to within 0.5 mag. 



6 THE RADIAL DISTRIBUTION OF 
SUBHALOS 

The large sample of subhalos in our simulations allows us 
to study their radial distribution and to investigate how 
it depends on the mass of the parent halo. In Fig. 6, we 
plot the number density of substructures as a function of 
the normalised distance R/R200 from the centre of the halo, 
defined here as the position of the most bound particle in 
the halo. We show averaged results for haloes with masses 
~ 10 15 h _1 M , ~ 10 14 h^MQ, and ~ 10 13 7i _1 M , and we 
limit the analysis to subhalos with masses greater than a 
fixed fraction (2 x 10 -4 ) of the virial mass of the parent 
halo. This fraction is chosen because it lies just above the 
mass limit where it is possible to identify substructures in all 
of our simulations. As Fig. 1 shows, there are typically ~ 50 
subhalos per parent halo with M sub /M2oo > 2 x 10~ 4 , so by 
stacking a large sample of haloes, it is possible to calculate 
an average density profile that has rather little noise. Note 
that the densities plotted in Fig. 6 have been normalised 
to the mean density inside the virial radius. The solid lines 
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Figure 5. Ratio in mass between the two most massive substruc- 
tures as a function of the parent halo mass. As in Fig. 4, small 
symbols represent the values measured for each individual halo, 
while the filled circles give the median in the same bins as in 
Fig. 4. The error bars mark the 20 and 80 per cent percentiles of 
the distribution. 



with symbols show results for the 3 different parent halo 
mass ranges defined above. For comparison, we have also 
plotted the dark matter radial profiles as dashed lines. Note 
that a small shift in the abscissa has been added to make 
the plot more readable. 

We find that the subhalo profiles are 'anti-biased' rel- 
ative to the dark matter in the inner regions of the haloes. 
This agrees with the results of Ghigna et al. (2000). Surpris- 
ingly we also find that the radial number density profiles are 
steeper in low mass haloes than in high mass haloes, a find- 
ing that reserves further investigation. 

We now use our highest resolution cluster simulation to 
investigate whether subhalos of different mass have different 
radial profiles. In Fig. 7, we show the cumulative fraction 
of substructures as a function of R/R200 for subhalos with 
M suh > 0.01 M200 (solid line) and M suh < 0.01 M 20 o (dashed 
line) . 

As Fig. 1 already made clear, there are many more 
substructures with M su b < 0.01 M200 than with M su b > 
0.01 M200 (9749 versus 96). Fig. 7 now shows that more mas- 
sive substructures are preferentially located in the external 
regions of their parent haloes. This can be understood as 
a consequence of tidal truncation and stripping effects that 
quickly decrease the mass of subhalos as they fall into the 
cluster and reach the dense inner cores of the parent haloes 
(see Section 7.2 for a more quantitative analysis of mass-loss 
due to stripping) . 

Also note that this finding can be naturally explained 
as a consequence of the orbital decay experienced by sub- 
structures. As shown in Tormen, Diaferio & Syer (1998), the 
orbital decay is consistent with expectations based on the 
combined effects of dynamical friction and mass-loss. As a 



Figure 6. Radial distribution of substructures in haloes of dif- 
ferent mass. Lengths are given in units of -R200 and densities 
are normalised to the mean density inside i?200- Symbols con- 
nected by solid lines show the number density profile of substruc- 
tures (filled circles are for 10 h~ Mq haloes, empty circles for 
10 14 /i _1 Mq haloes and filled squares for 10 13 /i _1 M haloes). 
Symbols connected by dashed lines show the corresponding dark 
matter profiles. 



result, massive substructures are driven to the centre more 
rapidly than less massive ones: Tormen et al. show that the 
orbital decay occurs in less than a Hubble time if the ini- 
tial mass of substructures is larger than 1 per cent of the 
mass of the main cluster, while the substructures can retain 
their identity for a significant fraction of the Hubble time 
if their mass is smaller than 5 per cent of the main clus- 
ter mass. Once driven to the centre, massive substructures 
are destroyed and no longer distinguishable from the central 
halo; this naturally explains the mass segregation that we 
see in our simulations. 

In Fig. 8, we show the cumulative fraction of the total 
mass of the parent halo that is in substructures as a function 
of normalised distance from the halo centre. The different 
lines and symbols have the same meaning as in Fig. 6 and 
represent median relations for all the haloes in each mass 
bin. The error bars mark the 20 and 80 per cent percentiles 
of the distribution. The mass fraction in substructures rises 
from ~ 1 per cent at a radius ~ 0.3 R200 to ~ 6 per cent at 
r ~ i?200- Note that the total mass fraction is dominated by 
the small number of most massive subhalos, and is hence a 
rather noisy quantity that shows large variations from sys- 
tem to system. The mass fraction may also reach values 
above ~ 10-15 per cent, but then the underlying FOF par- 
ent halo is typically quite aspherical, with the most massive 
subhalo often lying outside the formal radius R200 of the 
parent halo. 
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Figure 7. Cumulative radial distribution of the number of sub- 
structures in two different mass ranges in the simulation M3. The 
dashed line is for substructures with M su ^ < 0.01 M200, and the 
solid line is for M sub > 0.01 M200- 
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Figure 8. Cumulative mass fraction in substructures as a func- 
tion of the distance from the halo centre for haloes of different 
mass. The solid line is for haloes of mass ~ 10 15 /i -1 Mq , the 
dashed line for haloes of mass ~ lO 14 ft- 1 M , and the dash- 
dotted line for haloes of mass ~ 10 13 h _1 M0. The error bars on 
the last point mark the 20 and 80 percentiles of the distribution. 



7 SUBHALO HISTORIES 

So far we have analysed the properties of subhalos only at 
redshift 2 = 0. In this section, we turn to the time evolution 
of the masses of subhalos and their merging histories. 

In order to carry out this analysis, we have measured 
merger trees for each subhalo from the simulations. These 
trees allow us to specify all the progenitors (or the descen- 
dants) of a substructure at each epoch. To build the merger 
trees, we use a slightly modified version of the code described 
in Springel et al. (2001b). We briefly review the main fea- 
tures of the relevant algorithms in the following section. 



7.1 Constructing merging trees 

Following Springel et al. (2001b) we define a subhalo Sb at 
redshift z-g. to be the progenitor of a subhalo 5a at redshift 
za (with za < 2b) if more than one half of the Mink most 
bound particles belonging to Sb end up in Sa- Springel et al. 
(2001b) adopted iVn n k = 10. However, we obtained consider- 
ably better results with a value of JViink that varies between 
10 for the less massive substructures to 100 for more mas- 
sive ones. In this way, occasional failures of the code when 
building the merger trees were more robustly avoided. Par- 
ticularly if substructures undergo major mergers, the code 
identified in some cases the wrong progenitor for AT link = 10, 
or lost track of an entire subhalo. We have also updated 
the code so that volatile links to 'evanescent' substructures 
(i.e. objects close to the resolution limit that occasionally 
appear and then disappear) are avoided. 

With these choices, we manage to trace 85-90 per cent 
of all the substructures in our simulations back to the point 
when they were first accreted. This fraction goes up to 87-93 
per cent if we only consider substructures with more than 
100 particles. 



7.2 Mass accretion history 

We now use our merging trees to study the mass accretion 
histories of the subhalos in our simulations. Van den Bosch 
(2002) has carried out a similar analysis for dark matter 
haloes and has proposed an analytic expression for the mass 
accretion function based on the extended Press-Schechter 
formalism (Press & Schechter 1974; Bond et al. 1991; Bower 
1991). This function was found to be in excellent agreement 
with the results of high-resolution iV-body simulations. Our 
aim is to study the mass accretion function for subhalos 
and study whether there is any dependence on mass or on 
environment. 

We have selected subhalos at redshift 2 = in two 
different mass ranges (~ 10 11 /i _1 M© and ~ 10 12 h _1 M ). 
In order to test for the effects of environment, we selected 
on one hand subhalos that reside within the virial radius 
of the massive clusters that formed in simulations Bl and 
B2 by the present day, and on the other hand subhalos lo- 
cated within the smaller haloes found in simulation M3. In 
the following, we will refer to the substructures in the clus- 
ters as cluster subhalos and to the substructures inside the 
smaller haloes of M3 as field subhalos. Note that since we 
have excluded from our analysis the main subhalo associ- 
ated with the FOF group and since on average the most 
massive substructure has a mass a few per cent of M200 
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(see Fig. 4) we end up with very few substructures selected 
from M3. In particular we have only 5 substructures with 
a mass ~ 10 12 h~ 1 MQ and 38 substructures with a mass 
~ 10 11 ft~ 1 M Q . The corresponding numbers for the sub- 
structures selected from simulations Bl and B2 are 62 and 
338. 

For each of these samples, we build the mass accretion 
histories as follows: we start from a particular subhalo at 
redshift z = and construct its merger tree as described in 
the previous section. At each redshift we track the history 
of the selected subhalo by linking it with its most massive 
progenitor. 

In Fig. 9 we show a typical example of a mass accretion 
history for a subhalo with mass 2 x 10 11 /i _1 Mo at redshift 
zero. In the lower panel, we show the mass accretion history 
of the subhalo and in the upper panel, the corresponding 
mass of the halo in which the subhalo resides at each red- 
shift. 

In this example, the subhalo was accreted onto a larger 
halo at redshift ~ 1 (shown as a dotted line on the plot). 
For times prior to this event the substructure was a main 
subhalo, i.e. the subhalo corresponding to a FOF group, and 
its mass grew monotonically in time. From now on, we will 
refer to this event as the accretion time (t aC cr) of the subhalo. 
A few snapshots later, at redshift ~ 0.8, the substructure 
and its host halo were accreted onto the main progenitor of 
the cluster (shown as a solid line on the plot). 

After the subhalo is accreted, it suffers significant tidal 
stripping and decreases in mass. In this particular example, 
the final mass of the subhalo is ~ 40 per cent of the value 
at its accretion time. 

We find that for ~ 60 per cent of the subhalos in the 

10 11 mass bin and ~ 80 per cent of the subhalos in the 

10 12 bin the accretion event corresponds to the accretion 
of the substructure onto the main progenitor of the cluster 
itself. For most of the rest, the time elapsed between these 
two events is fairly short. The results we will show later 
are essentially unchanged if we adopt as definition of the 
accretion time, the accretion of the substructure onto the 
main progenitor of the cluster itself. 

Fig. 10 shows the distribution of the accretion redshifts 
for the cluster substructures in our sample. Interestingly, we 
find that a large fraction of the substructures are accreted 
at redshift z < 1. As noted above, for most of these sub- 
structures this accretion event corresponds to the accretion 
onto the main cluster itself. Our results hint that substruc- 
tures are constantly erased in the cluster, being replenished 
by newly infalling haloes. 

In Fig. 11, we plot the distribution of M(t = t )/M(t = 
£accr), i.e. the ratio of the mass of the subhalo at the present 
day to the mass it had when it was first accreted. The his- 
tograms show that this ratio has a quite broad distribution, 
varying from a value of ~ 1 to ~ 0.1. We note that most 
of the subhalos that have lost only small amounts of mass 
have been accreted very recently. 

This is more clearly shown in Fig. 12 where we plot 
the average mass accretion function for the cluster subhalos 
in the two mass bins considered. Three different accretion 
redshift intervals are considered and in all cases the sub- 
halo masses are normalised to the mass of the subhalo at 
£accr- The thick solid line shows the mean relation for sub- 
halos with mass ~ 10 11 h~ 1 M & , while the thin line shows the 
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Figure 9. Example for a typical mass accretion history for a sub- 
halo of mass 2 X 10 n /i Mq (lower panel), and the corresponding 
variation of mass for the parent halo in which the subhalo resides 
(top panel). The vertical solid line corresponds to the last time the 
subhalo is outside the main progenitor of the cluster; the dotted 
line corresponds to the time the subhalo becomes a substructure 
(see the text for details). 



relation for M aub ~ 10 12 ft _1 M . The mass accretion func- 
tion monotonically increases for times prior to the accretion 
event. Once the substructures are accreted, tidal stripping 
is effective and operates on short time scales. The longer the 
substructure spends in a more massive halo, the larger is the 
destructive effect of tidal stripping. Substructures remaining 
at z = that were accreted at redshift larger than 1 (panel 
c) have been typically stripped of ~ 80 per cent of their 
mass. There is also a slight indication that stripping is more 
effective for more massive substructures: panels (a) and (b) 
show that more massive substructures have been stripped 
significantly more than less massive substructures accreted 
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Figure 10. Distribution of the accretion redshifts for the cluster 
subhalos sample. A small shift is added to the abscissa to produce 
a more readable plot. 
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Figure 11. Distribution of the quantity M(t = to)/M(t = t aC cr) 
for the cluster subhalo sample. A small shift is added to the ab- 
scissa to produce a more readable plot. 



at the same redshift. This effect does not appear in panel 
(c) but note that we have very few objects accreted in this 
redshift bin in our sample of more massive subhalos. 

In Fig. 13 we compare the mass accretion histories of 
field and cluster subhalos. We limit the analysis to substruc- 
tures with mass ~ 10 11 /i _1 Mq. Again the mass accretion 
function is normalised to the mass of the subhalo at the 
accretion time. We find that field subhalos and cluster sub- 
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Figure 12. Average mass accretion history for 2 = cluster 
subhalos accreted in three different redshift bins. Thick lines are 
used for subhalos with mass ~ 10 11 (i^'Mq and thin lines are 
used for subhalos with mass ~ W 12 h- 1 M e . The histories arc 
normalised to the mass of the subhalo at the accretion time. 
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halos have remarkably similar histories suggesting that the 
efficiency of the tidal stripping is largely independent of the 
mass of the parent halo. 

7.3 Merging histories 

In hierarchical models of galaxy formation, galaxies and 
their associated dark matter haloes form hierarchically 
through merger and accretion processes. In this context, the 
term merger is usually used to refer to an interaction be- 
tween two objects of similar mass, while the term accretion 
is used to describe the infall of small objects onto much more 
massive haloes. 

Observational results and numerical simulations con- 
firm that interactions (such as tidal truncation or collisions) 
play an important role in the evolution of galaxies. There 
is, for example, solid evidence (Schweizer & Seitzer 1992; 
Whitmore et al. 1997; Barnes 1999) that at least some ellip- 
tical galaxies are the result of mergers between disk galaxies 
of similar mass. Mergers may also have a strong effect on 
the baryonic component of galaxies; they can trigger bar- 
instabilities in stellar disks and cause an inflow of gas into 
the galaxies centres, fuelling AGN activity or starbursts. 

Following the merger tree of the substructures in our 
sample we can analyse in detail the merger history, distin- 
guishing between mergers and accretion events. To build the 
merger history we proceed as follows: we consider all the sub- 
structures within the virial radius at redshift zero and follow 
their merger trees back in time, checking as substructures 
are accreted onto the main progenitor. We count as mergers 
all accretion events involving haloes with mass larger than 
2 x 10 10 tr 1 Mq and mass ratio smaller than 5:1. Note 
that the lower limit on the mass of the merging haloes cor- 
respond to the resolution limit of our cluster simulations (see 
Table 1). 

In order to have enough information without running 
into numerical resolution effects, we limit the present analy- 
sis to subhalos with mass ~ 10 12 h^ 1 M Q . Fig. 14 shows the 
mean number of mergers per subhalo, identified at redshift 
z = 0, occurring after the redshift plotted in the abscissa. 
As in the previous section, the sample is split into three 
different accretion redshift intervals. 

Merger events are less frequent once a halo is accreted 
onto a more massive structure. This is because the merging 
efficiency is higher in environments where the relative veloc- 
ities of subhalos are similar to their internal virial velocities. 
Once haloes are accreted by a massive halo, merging is sup- 
pressed by the large velocity dispersion they aquire. This 
effect is clearly visible in Fig. 14 in the change in slope of 
the curves near the accretion redshift. Note however that 
when one goes to significantly higher redshift the difference 
between the three curves vanishes. 

In Fig. 15 we show the fraction of subhalos that have 
had at least one merger after the redshift plotted in the 
abscissa. Note that the merger events we are considering 
are characterised by similar masses and will most likely in- 
fluence the morphology of the main galaxy, leading to the 
formation of a bulge component. The final morphology of 
the galaxies residing in these substructures will depend on 
the time between the last major merger and accretion onto 
the cluster: the longer this time, the larger is the likelihood 
that the galaxy can grow a new disk. 
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Figure 13. Mass accretion histories for our lowest mass bin 
(~ 10 11 /j^'Mq) for cluster subhalos (solid line) and for field 
subhalos (dashed line). As in Fig. 12, the sample is divided into 
three different subsamples according to the redshift of accretion. 
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Figure 14. The mean number of mergers after rcdshift z for 
substructures selected from the simulations Bl and B2 at red- 
shift z = and with mass ~ 10 12 h~ 1 M.Q. The solid line is for 
substructures accreted at rcdshift z < 0.5, the dashed line for 
substructures accreted at rcdshift 0.5 < z < 1.0 and the dashed- 
dotted line for substructures accreted at rcdshift z > 1.0. 

The results in Fig. 15 show that ~ 80 per cent of the 
~ 10 12 h -1 M© haloes in our 2 = sample have had at 
least one major merger at redshift below 2; this fraction 
decreases to ~ 40 per cent for redshift < 1. Surprisingly the 
fraction is almost independent of the accretion time. These 
results suggest that a large fraction of subhalos in our cluster 
sample will host early-type galaxies with a significant bulge 
component. 



8 SUMMARY AND DISCUSSION 

We have used a set of high resolution numerical simulations 
in a ACDM Universe to study dark matter halo substruc- 
tures. Such dark matter substructures mark the sites where 
luminous galaxies are expected to be found, so the anal- 
ysis of their mass functions, radial distributions, merging 
and mass accretion histories should help us to better un- 
derstand the properties of the galaxies that form in hierar- 
chical galaxy formation models. Comparison with observa- 
tional data should then suggest the physics that needs to be 
included in viable models of galaxy formation and evolution. 

In agreement with previous work (Moore et al. 1999b) , 
we find that the shape of the subhalo mass function is almost 
independent of the mass of the parent halo, with galactic 
haloes being essentially scaled versions of cluster haloes. 

We find that the average mass of the largest substruc- 
ture within the virial radius (excluding the BCG) scales lin- 
early with the mass of the parent halo. If the stellar masses 
scale linearly with the dark matter mass of the parent sub- 
structure, the second ranked galaxies should have K-band 
luminosities that increase roughly linearly with the mass of 
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Figure 15. Fraction of substructures with at least one merger 
event at redshift z m < z. Different line styles are for different 
accretion redshift as in Fig. 14. 

the main halo and are equal to those of the third ranked 
galaxies to within 0.5 mag. 

Note that the assumption that the stellar mass scales 
linearly with the dark matter mass of the parent substruc- 
ture cannot generally be true since stars are typically much 
more concentrated than dark matter. The relation between 
stellar mass and the mass of parent substructure is then 
quite complex and should be followed considering the de- 
tails of star formation and feedback process as is done for 
example in Springel et al. (2001b). 

We have also used the simulations to study the radial 
distribution of substructures. In agreement with previous 
work (Ghigna et al. 2000), we find that the radial profile 
of substructure number density is 'anti-biased' relative to 
the dark matter profile in the inner regions of haloes. The 
most massive substructures reside preferentially in the outer 
regions of haloes. This is, at least in part, because substruc- 
tures undergo substantial tidal stripping in the dense inner 
regions of haloes. 

We have studied the evolution with time of this strip- 
ping process and find that once a subhalo is accreted by a 
larger system, tidal stripping is highly effective; the longer 
a substructure spends in a more massive halo, the larger is 
the destructive effect. This suggests that substructures are 
constantly erased in clusters, being replenished by newly in- 
falling galaxies. 

Interestingly, a significant fraction of the substructures 
found in present-day clusters were first accreted at redshifts 
z < 1, implying that tidal truncation of the dark haloes 
of cluster galaxies happened relatively recently, and that at 
z > 1 many present-day 'passive' cluster members should 
still have been central galaxies of their own extended dark 
haloes. If gas was able to cool in these haloes, the galaxies 
may have been considerably more active at z ~ 1 than at 
present. Substructures in smaller haloes have histories re- 
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markably similar to those in cluster haloes, suggesting that 
the efficiency of tidal stripping is largely independent of the 
mass of the main halo. 

Once a substructure is accreted onto a cluster, its merg- 
ing probability decreases because of the large velocity dis- 
persion of the system. Observational data are currently be- 
ing collected on merger rates in different environments (van 
Dokkum et al. 1999; Patton et al. 1997). In future work, 
we plan to carry out a more direct comparison with obser- 
vational data of this kind, and so to test the hierarchical 
paradigm for galaxy formation in new ways. 
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